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Optimal feedback control of turbulent channel flow 

By Thomas Bewley, Haecheon Choi, Roger Temam 1 , AND Parviz Moin 

Feedback control equations have been developed and tested for computing wall- 
normal control velocities to control turbulent flow in a channel with the objective of 
reducing drag. The technique used is the minimization of a “cost functional” which 
is constructed to represent some balance of the drag integrated over the wall and 
the net control effort. A distribution of wall velocities is found which minimizes this 
cost functional some time shortly in the future based on current observations of the 
flow near the wall. Preliminary direct numerical simulations of the scheme applied 
to turbulent channel flow indicates it provides approximately 17% drag reduction. 
The mechanism apparent when the scheme is applied to a simplified flow situation 
is also discussed. 

1. Motivation and objectives 

It is the goal of this project to study methods to counteract near-wall vortical 
structures in turbulent boundary layer flow using an active control system in an 
effort to reduce drag. From this study, we hope to better understand the physics of 
drag producing events and the sensitivity of boundary layer flow to control. As a 
more far-reaching goal, we would like to better understand how to develop control 
equations for general flow control problems, utilizing the equations governing fluid 
flow to achieve performance that is in some sense optimal for a given situation. 

With a well-chosen scheme using wall control only, it has been shown that a 
turbulent flow may be smoothed out in a near-wall region, and the drag may be 
substantially reduced. This scheme applies small amounts of wall-normal blowing 
and suction through the computational equivalent of holes drilled in the wall. Pre- 
vious ad hoc schemes by Choi et al. (1992) have reduced the drag by as much as 
20% by countering the vertical velocity slightly above the wall with an equal but 
opposite control velocity at the wall. The objective of this work is to derive more 
effective schemes by applying optimal control theory, utilizing the equations of mo- 
tion of the fluid to reveal the dominant physics of the control problem and the most 
efficient distribution of the control energy. This work is an outgrowth of the work 
done by Choi et al. (1993), where optimal control theory was applied to the stochas- 
tic Burgers equation. Here, we apply the theory to the Navier-Stokes equations, 
which necessitates a more involved treatment of the equations and more extensive 
computer resources. The scheme discussed in this report depends on measurements 
of flow velocities above the wall — this is not feasible in a practical implementa- 
tion. The scheme will later be reduced to a more practical one involving only flow 
quantities which are most easily measured in an experimental rig. 
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The model problem we study in this work is the turbulent flow inside a small 
segment of a fully developed turbulent channel (i.e. flow between two parallel walls, 
far from the inlet). This flow is governed by the same vortical structures as turbulent 
boundary layer flow in the near- wall region. 

Thus, the problem under consideration is a turbulent channel flow with no-slip 
walls and wall-normal control velocities <f > . Control will be applied to this flow in 
order to decrease the drag integrated over the walls at the expense of some measure 
of the net control effort. A feedback control algorithm has been developed which 
minimizes a “cost functional” constructed to represent this balance of the drag and 
the control effort. This method is introduced in Section 2. The control equations 
have been coded and tested in a direct numerical simulation of turbulent channel 
flow. Section 3 discusses preliminary results of these calculations. 


2. Formulation 


2.1 State equation (Navier- Stokes equation) 

As described above, the problem under consideration is a constant-flux turbulent 
channel flow with no-slip walls and wall-normal control velocities <j>. This problem 
is governed by the unsteady, incompressible Navier-Stokes equation, the continuity 
equation, and a constant flux integral constraint equation inside the domain ft and 
appropriate boundary conditions on the walls w (periodic conditions are implied on 
the remainder of the boundary of the domain T): 


duj d _ dp 1 d d 
&t + dxj UjUl dxi + Re dxj dxj U 1 



ix i dx\ dx 2 dx 3 


C 


in f 2 


(la) 

(lb) 

(1 c) 


uj = 0 j 

U2 = <f> > on walls, 

U3 = 0 j 


( 2 ) 


where x\ is the streamwise direction, x 2 is the wall-normal direction, .r 3 is the 
spanwise direction, Ui are the corresponding velocities, and p is the pressure. The 
constants in the problem are C (a measure of the flux in the channel) and Re (the 
Reynolds number). 


2.2 Optimal control of state equation 

The goal of controlling the channel flow is to minimize the drag on a section of 
wall with area A over a period of time T using the least amount of control effort 
possible. The relevant quantities of interest are thus the time averaged drag 



du\ 

dn 


dx x dx 3 dt 


( 3 ) 
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(where n is a unit vector in the inward wall normal direction) and a term represent- 
ing the expense of the control. The latter term may be taken to be the integral of 
the magnitude of the power input, which may be written 

T 

J JJ 1^ + P^ 2 / 2 )| dxi dxz dt, (4a) 

In addition, depending on the physical mechanism used to provide the control veloc- 
ities, the rate of change of the control hardware settings might be another important 
expense (for instance, representing the expense involve in changing the settings of 
control valves in the system): 

(46) 

A physically appropriate cost functional for this problem, then, balances the expense 
of the input versus the drag: 


•MJL\ 


d4> 

di 


dx i dx 3 dt. 


J{4>) — E\ + ^2^2 + D , ( 5 ) 

where t\ and I 2 are appropriate weighting factors. We could proceed from this 
point to attempt to construct a control procedure designed to minimize this cost 
functional. A mathematically more simple cost functional for the purpose of control 
theory (for reasons which will become evident as the control equations are derived) 
is quadratic in <f>. Physically, this represents the integral of the magnitude of the 
kinetic energy per unit mass input to the system, and may be written 

m = l Jr 1 JJj’ dx “ ,x,d, + ifl jj m l£ dx ' dx ‘ dt - (6) 

It will be seen later that, in most problems that we consider, the expense terms are 
much less significant than the drag terms (in other words, the control is relatively 
cheap). The use of other expense terms does not cause much additional complexity 
or insight into the method. 

The optimal control procedure considered, then, involves reducing the cost func- 
tional (6) for some period of time T. This method is described in Abergel and 
Temam (1990) in a related situation and is also discussed in Lions (1969). How- 
ever, this is a prohibitively expensive procedure for present computational resources 
because it involves storage and manipulation of several three-dimensional fields over 
the entire time period under consideration. The complexity of such an algorithm is 
discussed further in Choi et ai (1993). 

We therefore resort to a suboptimal control procedure (Choi et al 1993). In this 
method, the state equation is discretized in time, then a control procedure is applied 
to reduce an instantaneous version of the cost functional (6) 

= 2l Ji *' dx ‘ dXx + \JLl£ dx ' <ttI < 7 > 
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at each time step. 

By applying the control at each time step, the algorithm gives the control which 
minimizes the cost functional over some short time interval. Note, however, that 
this method does not look ahead to anticipate further development of the flow, and 
thus the solution by this method does not necessarily correspond to the solution by 
the optimal control method. Thus, posing the problem in this suboptimal form is 
another level of approximation to the physical problem of interest. 

The differences in complexity between the optimal and suboptimal schemes de- 
scribed above may be realized by drawing an analogy to a computer algorithm to 
play chess. A suboptimal chess program looks ahead one step to determine the 
move that leaves as good a position on the board as possible. Similarly, a subop- 
timal turbulence control scheme looks ahead one time step to determine the set 
of control velocities that leaves as good {i.e. low) a value of the cost functional as 
possible at the next time step. An optimal chess program, on the other hand, in- 
vestigates all possible developments of the game a certain number of steps into the 
future (knowing how the other player may respond), and then moves in the direction 
that leads to the best final position on the board. Similarly, an optimal turbulence 
control scheme investigates all possible developments of the flow a certain amount 
of time into the future (knowing approximately how the flow will respond), and 
then applies the set of control velocities that leads to the best (i.e. lowest) time- 
averaged cost functional. Such a method requires significantly more resources than 
the suboptimal method. 

2.8 Time discretization of state equation 

The suboptimal control procedure introduced above is now applied to the state 
equation (1). To do this, we discretize (1) in time, then apply a feedback control 
algorithm to modify the flow at the next time step. A consistent approach is to 
use a second order Crank- Nicolson method (implicit) on all terms. The momentum 
equation (la) thus takes the form: 


u: — u 


n — 1 


At 




"- 1 U «- 1 ) = 


l(dp n dp^_ \ 1 ( 

2 \ dxi dx, ) 2Re \ 


d d 


< + 


d d 


dxj dx } * dx j dx 




( 8 ) 


where a superscript n indicates the value at time step n. 

It is now useful to put the time discretized form of the entire state equation 
governing the flow in the domain into the form 

£ n + K n - 1 =0, (9a) 


where K n contains all the terms which in some way depend on the state variables 
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from the current time step, and 7Z n 1 contains the remaining terms: 
d d / dp' n dP " d \ 

<-e'er i er i < + K-k + ir, s ‘' + ar^<) '* 




in (96) 


ft 


UL« 


dx i dx 2 dx 3 


* ( n— 1 q & & n-1 

- U t - Pi u. 


n"- 1 = i 


dxj dxj ' 


W^+r'^r'-r 1 ) 


0 

C' . 


in Q 


in Q 


(9c) 


In the above equation, j3\ = At/ 2 Re, fa = At/2, dP/dx \ is the mean pressure 
gradient in the x\ direction (adjusted at each time step to provide constant mass 
flux), and p ' accounts for the pressure variations within the domain (periodic in x\ 
and 23 ). Note that (lb) and ( 1 c) have been multiplied by constants to obtain ( 9 ). 

Associated with this problem are the boundary conditions B: 

B\ = u\ = 0 

I?2 = U 2 = <j> (10) 

Bj — U 3 — 0 . 


The “flow problem”, which will hereafter be denoted £/ , is taken to refer to the 
differential equation (9) together with the boundary conditions (10). 

2.4 Suboptimal control of state equation 

In this section and the next, we develop a method to solve for the gradient of 
the cost functional J and with this a control procedure based on this gradient 
information to minimize J at each time step. 

Consider the Frechet differential (Vainberg, 1964) of the cost functional J in (7): 

_ lim J(<t> + 4)- J(<t>) 

= 1 Sfj idl ' *' 3 + 2 JI. iL (W*) iX ' iX> ' 


The gradient of the functional J with respect to the control distribution <j> may be 
extracted from this equation by expressing the last term on the RHS in terms of an 
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inner product on <j>. It is for this reason that we now formulate what we shall call 
the “differential problem”. 

Define 0 using a Frechet differential such that 


= = l im 

~D 4 > «— • -o e 


( 12 ) 


where <f> is some arbitrary or “test” distribution of control velocities. Thus, 0" 
is a differential state representing the sensitivity of the state U n to control for a 
particular control distribution <j> n applied over the time duration (< r<_1 , t n ). The 
differential 0 is decomposed into components in a fashion similar to the state 

( Ui ( ri , x 2 , r 3 )\ / 9,{x 1 ,x 2 ,x 3 ) 

p'(x i,x 2 ,x 3 ) J, 0=1 p(xi,x 2 ,x 3 ) 
dP/dx j / \ A 

The equations governing the differential state 0” follow directly by taking the 
Frechet differential of the state equation (9) and its boundary conditions (10). Note 
that the term TZ n ~ l in (9) does not depend on 4> n and thus makes no contribution. 
The contribution from the term AC" is linear and may be written 



A"0" = 0, 


(13a) 


where 


AO = < 


* _ A + A (^ + ““ + + Ui w) 


dp 


dui dOi 


ddj 

dxj 


W//„ 


$i dx\ dx 2 dx 3. 


The boundary conditions B, from (10), are 


b x = e x = o 
J3 2 = — i 

— $3 = 0 . 


in fl 
in Q 


(136) 


(14) 


The differential problem”, which will hereafter be denoted ^/, is taken to refer to 
the differential equation (13) together with the boundary conditions (14). 

Consider again the Frechet differential of the cost functional J in (11): 

= 2 JJjid x. ^ + 2 II §£ *■. *»■ (15) 

The gradient of the functional J with respect to the control distribution <p may be 
extracted from this equation by expressing the integral of 06 \ jdn in terms of an 
inner product on <£. This may be done by solving the differential problem «c/, as is 
done below. 
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2.5 Solution of differential problem &/ by adjoint method 

An “adjoint problem” is now formulated which may be used to bypass direct 
solution of the differential problem £/ itself. 

Define an adjoint operator A* using the equation 


< A0, V > = < 0, > + 6, (16) 

where A (which depends on U) is defined in equation (13b), the boundary conditions 
on 0 are given in equation (14), and an adjoint state $ has been defined in a fashion 
similar to U and 0: 

( ^•(xi,X2 ,x 3 )\ 
n(xi,x 2 ,X3) I . 

The adjoint operator is formed by moving all of the derivatives in the inner product 
(the integral over the volume of the product of the two terms, denoted <*,•>) 
from the differential 0 to the adjoint It is a straightforward exercise to write 
out the volume integrals corresponding to the LHS of (16) and then to rearrange 
this expression into the form of integrals corresponding to the RHS of (16) using 
integration by parts. From this is deduced A* and the condition at the boundary 
resulting from the wall terms, which are all placed into the expression for b : 

b =< A0,$ > - < 0,A*$ > . (17) 

Through equation (13a), the first term on the RHS of equation (17) is zero. If we 
form a similar homogeneous adjoint differential equation for the adjoint $ 


A** = 0, 


(18) 


with boundary conditions as yet undetermined, then equation (17) reduces to 


6 = 0 . 

Using the method described above, it is easy to show that 

d d , . „ / < 97 t „ duj drpi 


f *" A ^a^* +A (s; + ‘‘ <i, + *>5i 


A*tf = { 


(sX j 


in ft 


and 


[ 02 JJJ n \j>\ dx\ dx 2 dx 3 
-//.(*£< 

+ 4>(p2 »2 » ~ fll - $2 <t>K2 ^ 2 )^ dx j dx 3 = 0. 


89 

•(V’l) - h pni{V>2) + A^Ofo) 


(19) 


(20) 


( 21 ) 
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(Note by comparison of (20) with (13b) that the operator A is not self-adjoint due 
to the effect of the convective terms of the momentum equation.) These adjoint 
equations may be exploited to solve the differential problem srf . 

We now formulate an “adjoint problem”, which will hereafter be denoted 
defining an adjoint state with the homogeneous differential equation (18) and 
with accompanying boundary conditions B* as yet undefined. Note by the above 
discussion that one of the by-products of the formulation of this problem is the 
relation at the boundary given by (21). We are now at liberty to choose boundary 
conditions for the adjoint problem such that this relation is useful — it is exactly 
for this reason that the formulation of an adjoint problem is considered. With this 
in mind, we may choose the boundary conditions B* as 


B\ = 0i = l 

B * = 0 2 = 0 ( 22 ) 

B ; = 03 = 0. 


Using these boundary conditions and the continuity equation for the adjoint velocity, 
equation (21) reduces to 



d6 x 

dn 


dx i dx 3 


= 2 


JL 


<f> Ren dx\ dx 3. 


(23) 


To compute the RHS, we must solve the adjoint problem £f *. This is done numeri- 
cally and must be repeated at each time step as A* changes as the flow U develops 
with time. 

The differential of the cost functional (11) may be rewritten using (23) as 


^ JJ <t> 4> dxi dx 3 - n ? - ~ - JJ nfidxi dx 3 , (24) 

where 7r is the adjoint pressure on the wall. Finally, the desired gradient of the cost 
functional J may be extracted (Vainberg, 1964): 


VJ{<p) £ U 2 Re 

— — 0 — 7T. 

V<f> A v A 


(25) 


A feedback control procedure using a simple gradient algorithm at each time step 
may now be proposed such that 


0 n ’* +1 - <j> n * k 


VJ(cf> n ’ k ) 

p V4> 


(26) 


where superscript n indicates the time step as before and superscript k indicates 
an iteration step at that particular time step. This algorithm attempts to update 
0 in the direction opposite to the local direction of increase of J . For fixed n as 
k — ► oo with sufficiently small /j, this gradient algorithm should converge to some 
local minimum of J over the control space 0 if the approximation of VJ /P0 is 
sufficiently accurate. Note, however, that as the time step n advances, J will not 
necessarily decrease (Choi et al, 1993). 
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3. Accomplishments and future work 

S.l Elementary drag reducing mechanisms 

Choi et al. (1992) found that by applying a control velocity equal and opposite to 
the vertical velocity at = 10, a drag reduction of nearly 20% could be achieved. 
Vertical transport of streamwise momentum in the near-wall region (primarily due 
to longitudinal vorticity) produces “sweep” events and thus local regions of very 
high drag. Applying a countering control velocity tends to reduce this effect. A 
related mechanism described by Lumley (1993) further explains these results; con- 
trol applied to reduce the spinning of the near-wall vortices reduces their energy, 
stabilizing them in space and thereby reducing the “bursting” frequency, which also 
tends to reduce the drag. 

In the tansverse plane, countering the vertical velocity above the wall corresponds 
to a control which de-spins the near- wall vortices, as shown in Figure 1. This process 
leads to the removal of fluctuations in the near-wall region, which diminishes the 
mixing capability of the turbulence and therefore reduces drag. This type of control 
corresponds to blowing where the drag is high, which decreases the high velocity 
gradients at the wall and thus smooths out the flow in the near- wall region, as 
shown in Figure 2. 



FIGURE 1. Stabilization mechanism in cross flow plane. The effect of the control 
velocities shown is to de-spin the near-wall vortex, reducing momentum transport 
near the wall. 



FIGURE 2. High drag is decreased by blowing at the expense of suction in the 
regions of low drag, resulting in a net smoothing of the near-wall velocity profiles. 
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Figure 3 shows the application of the suboptimal control scheme to a simple flow 
configuration of longitudinal vortices embedded in an initially parabolic flow. A 
cross flow plane is shown. In regions below downward moving fluid (sweep events) 
the streamwise (into the page) drag is higher and blowing is applied. In regions 
below upward moving fluid (ejection events), the streamwise drag is lower and 
suction is applied. The overall control distribution from the suboptimal scheme is 
in a sense that acts to de-spin the near-wall vorticity, and thus acts in accordance 
with mechanisms described above. 



u = 0 


Control velocities 


FIGURE 3. Optimal control scheme applied to longitudinal vortices. Interior 
vectors are cross flow velocities and contours are of streamwise velocity, indicating 
a sweep event between two near-wall vortices and ejection events outside of them. 
Control velocities shown on the wail (not to scale) indicates blowing at the sweep 
event and suction at the ejection events. 

The adjoint analysis utilizes all the information present in the near-wall region to 
extract the sensitivity of the instantaneous drag to the variation of the control. This 
scheme may be reduced to an approximate one relying only on wall information by 
approximating the near-wall velocities using a Taylor’s series extrapolation of the 
velocity gradients at the wall. The correlation between the full adjoint analysis and 
approximations of the adjoint problem using only information available at the wall 
is still being investigated; preliminary results indicate that the performance is not 
severely degraded by this approximation. 
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3.2 Suboptimal control of turbulent channel flow 
The scheme introduced in Section 2 was tested by applying it to a direct numerical 
simulation of turbulent channel flow. A 17% drag reduction was seen as compared 
to a flow with no control. Results are plotted in Figure 4. This calculation was 
done in a flow with Re r — 100 based on the friction velocity and the channel half 
width using a 32x65x32 grid and the spectral method of Kim et al. 1987. Although 
these results should be considered preliminary, they are quite promising. 



I 1 i r i 1 

0 10 20 30 40 50 


time 

FIGURE 4. Performance of suboptimal scheme compared to no control and the 
scheme of Choi et al. (1992). Parameters for suboptimal scheme are fi = 0.01, 

t = 10, T + — 1. Legend: suboptimal scheme, <j> = —^+= 10 ^ 

no control. 


3.3 Future work 

At present, the drag reduction obtained using a suboptimal control scheme is 
still slightly less than the drag reduction obtained using the ad hoc scheme of Choi 
et al. (1992), as shown in Figure 4. It is hoped that by further variation of the 
parameters and careful study of the numerical issues of the adjoint pxoblem, the 
result using the suboptimal formulation may be significantly improved. We expect 
that, using the suboptimal method, a significant improvement is possible over all 
ad hoc schemes, as the suboptimal scheme uses the entire flow information in the 
near- wall region and is rigorously based. Also, work is currently in progress with 
Dr. Chris Hill to reduce the suboptimal control scheme to one which depends on 
wall information only. Preliminary results of this work are also quite promising — 
a discussion of this project is included in the next report in this volume. 
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